A multiscale study of electronic structure and quantum transport in C 6n 2 _£f 6n -based 

graphene quantum dots 
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Wc implement a bottom-up multiscale approach for the modeling of defect localization in C 6n 2 Hen 
islands, i.e. graphene quantum dots with a hexagonal symmetry, by means of density functional 
and semiempirical approaches. Using the ab initio calculations as a reference, we recognize the 
theoretical framework under which semiempirical methods describe adequately the electronic struc- 
ture of the studied systems and thereon proceed to the calculation of quantum transport within 
the non-equilibrium Green's function formalism. The computational data reveal an impurity-like 
behavior of vacancies in these clusters and evidence the role of parameterization even within the 
same semiempirical context. In terms of conduction, failure to capture the proper chemical aspects 
in the presence of generic local alterations of the ideal atomic structure results in an improper de- 
scription of the transport features. As an example, we show wavefunction localization phenomena 
induced by the presence of vacancies and discuss the importance of their modeling for the conduction 
characteristics of the studied structures. 



I. INTRODUCTION 



Graphene is a carbon allotrope material that has trig- 
gered a vast interest within both academic and indus- 
trial communities for its peculiar electrical, mechanical 
and optical characteristics [l|. Additionally, graphene re- 
sults extremely convenient in terms of electronic struc- 
ture modeling, since its planar monolayer topology and 
sp 2 hybridization allow for a simple and accurate ir- 
orbital tight-binding (TB) description of its electronic 
bands @, y|. This approach has been extensively imple- 
mented also in the case of quantum transport studies of 
quasi one-dimensional graphene constrictions with arm- 
chair or zigzag edges, widely known as graphene nanorib- 
bons (GNRs)[4j-l7|. In the course of these few years since 
graphene's laboratory isolation and due to the intense 
research on the field, experimental knowledge has grown; 
e.g. recent transmission electron microscopy images of 
monoatomic graphene flakes obtained with mechanical 
exfoliation have shown Stone- Wales defects and vacan- 
cies on the crystal membrane Q, while chemical reac- 
tivity with metal oxides has been achieved both in de- 
fected and in edge sites 0. Theory on the other hand by 
means of first principles calculations has demonstrated 
energetically favorable chemisorption processes on the ar- 
eas that diverge from the perfect two-dimensional atomic 
lattice [Tot- The aforementioned advances, among oth- 
ers, have raised new challenges in the field of graphene- 
based modeling in chemically and geometrically com- 
plex environments. Ideally, density functional theory 
(DFT) with appropriate exchange-correlation function- 
al could formulate an accurate and transferable framc- 
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work for simulations on these systems [Til , [ill [l2j . The 
undoubted intrinsic rigor of ab initio approaches comes 
in hand though with computational overload limitations. 
In this sense, only order ~ 10 3 can be affordable for all- 
electron bandstructure calculations (under linear scaling 
optimization techniques fl3l - [l6l ] ) . while an extra load has 
to be considered for self-consistent nonequilibrium quan- 
tum transport studies. Therefore, it becomes evident 
that a massive atomic reconstruction study in the pres- 
ence of generic local alterations of the symmetry or the 
ideal atomic structure cannot be addressed with 'dog- 
matic' singular theoretical approaches (e.g. pure ab ini- 
tio methods, pure TB studies etc.), whereas a versatile 
multilevel approach would seem more appropriate. 

In this article a multiscale study of electronic struc- 
ture and quantum transport is carried out for graphene 
quantum dot complexes. The study's groundwork fo- 
cuses on: (i) the structural characteristics, and (ii) the 
methodological approach. The islands under consider- 
ation are the coronene molecules [13, EH with a general 
chemical type of C§ n iH§ n in a pure, defected (with a sin- 
gle vacancy) and hydrogen functionalized form. These 
can be thought of as planar complexes of benzene rings 
that grow rotating around a central benzene ring, form- 
ing six hydrogen-passivated zigzag edges (see fig. [T]). 
Methodologically, electronic structure is initially studied 
with DFT while optical properties are calculated within 
the time-dependent density functional theory (TD-DFT). 
These results serve as a reference for the comparison 
with similar calculations by means of two parameter- 
ized semiempirical methods: (i) the extended Hiickel 
(EH) theory [l9| and (ii) , the next- neighbor tight-binding 
model. As soon as a proper functional framework is iden- 
tified for the semiempirical approaches the study pro- 
ceeds with the calculation of electronic transport. Par- 
ticular attention is paid to the effect of defect localiza- 



2 




FIG. 1: C 6n 2H(,„ molecular complexes: a) C24-H12 (coronene, 
n=2), b) C54H1S (coronene 19, n=3), c) C96-H24 (coronene 
37, n=4), d) C150-H30 (coronene 61, n=5) 



tion both within a level of characterization as well as 
model calibration, since the computational results indi- 
cate impurity-like behavior of single vacancies in these 
systems. Conceptually, although the final objective is 
physical (quantum transport modeling in graphene is- 
lands), the basis is founded on chemistry (multiscale com- 
parative analysis of the electronic structure). 

The paper is organized as follows: In section [TTI we re- 
view the methodological steps of the theoretical model, 
section IIIII presents first principles and semiempirical 
electronic structure calculations for C§ n iH§ n islands up 
to n = 5, section IIVI analyzes transport characteristics 
and focuses on the effects of local atomic reconstruction 
on the conductance, while in section [V] we discuss our 
results. 



II. METHODOLOGY 

Geometry relaxation and electronic structure proper- 
ties (eigenvalues, eigenfunctions, density of states) of 
various C§ n iH§ n clusters are extrapolated by DFT cal- 
culations on a split- valence double-zeta (3-21g (g^-HH) 
and a minimal (STO-3G [HI, H3) basis set, as imple- 
mented in the GAUSSIAN code[25|]. The semiempirical 
three-parameter hybrid nonlocal exchange and correla- 
tion functional of Becke and Lee, Yang and Parr [2^- 
f29( | (B3LYP) has been chosen here for its capacity to 
predict a large range of molecular properties for aro- 
matic systems [30l 131] . Additional optical properties (ex- 
citation energies, fundamental optical gaps) are calcu- 
lated within a TD-DFT approach for comparison be- 
tween theory and experiment. Electronic structure re- 
sults are then confronted with similar ones obtained by 
two semiempirical methods [H, (HI that also present a 
precision/efficiency mismatch among them: (i) the ex- 
tended Hiickel method, and (ii) a next-neighbor tight- 
binding model. In the case of the EH method three dis- 



tinct parameterizations are used: (i) the first one con- 
siders a standard valence 2s2p-basis set of single-^ Slater 
orbitals for C atoms, principally deriving from the initial 
values used by Hoffmann (l9j (EH-sp from now on)[48j. 
(ii) The second one is a 2s2p3d-based parameterization 
with valence/polarization double-C exponents and C pa- 
rameters fitted to recreate the bandstructure of two- 
dimensional graphene as given by DFT calculations [33- 
36] (EH2-spd from now on) . (iii) The third parameteriza- 
tion derives in the similar way to the second o ne, whereas 
here the polarization orbitals are absent [3J, [35|(EH2-sp 
from now on). For the first- neighbor TB model a stan- 
dard to = 2.7eV C-C hopping integral is used, while va- 
cancies are approximated with the insertion of a local 
point potential U — > 00 (unless explicitly referred to in 
the text) . Although all methods construct the molecular 
orbitals on the basis of the linear combination of atomic 
orbitals there is a distinct difference in the level of accu- 
racy that each method delivers. In the DFT case the ba- 
sis set is comprised of Gaussian-type orbitals with weight- 
ing coefficients that are both calculated self-consistently 
in order to reproduce the best approximation of the exact 
ground state density of the system. In the EH case the 
bases are nonorthogonal Slater-type orbitals with fixed 
weighting coefficients that have been parameterized on 
the basis of experimental data or first principles calcula- 
tions. Finally in the TB case no real orbitals exist and 
the system Hamiltonian is constructed by a simple finite 
difference approach on next-neighbor atoms that through 
a proper choice of the hopping integral is representative 
of the 7r-orbital in the sp 2 hybridization scheme. Natu- 
rally, the level of computational efficiency is the inverse, 
ranging from molecular (DFT) to mesoscopic (TB). 

Electronic structure results are obtained through a di- 
rect diagonalization of the respective Hamiltonian ma- 
trix. Comparisons take place in terms of highest occu- 
pied molecular orbital (HOMO) and lowest unoccupied 
molecular orbital (LUMO) gaps, energy eigenstates e a 
and their respective eigenfunctions ty a . The local den- 
sity of states LDOS(r, E) at the positions r of the device 
atoms at energy E is calculated as: 

LDOS(r,E) = Y,\* a (r)\ 2 S(E-e a ), (I) 

a 

where 8 is the Delta function, while summing over all 
atoms gives the total density of states (DOS) of the 
molecular systems at this energy. 

Quantum transport is calculated within the non- 
equilibrium Green's function formalism [37j. In partic- 
ular, the method is based on the single particle retarded 
Green's function matrix G = [ES — H — — 
where E is the energy, H and S are the device Hamil- 
tonian and the overlap matrix respectively (written in 
an appropriate basis set), while T^l,r are the self-energy 
matrices that account for the effect of scattering due to 
the left (L) and right (R) contacts. In the TB case 
the overlap matrix coincides with the unitary one. The 
^l,r terms can be expressed as S = Tg s r\ where g s is 
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the surface Green function specific to the contact type 
and r is the Hamiltonian relative to the interaction be- 
tween the device and the contact. The calculation of 
the Green's function permits for the evaluation of all 
the quantities of interest for conduction, e.g. the de- 
vice spectral function is the anti-hermitian part of the 
Green's function A — i(G — G^) from which the Den- 
sity of States can be obtained as D(E) = ^Trace(AS). 
Moreover in the coherent transport regime, the expres- 
sion used for the zero-bias transmission probability reads 
T{E) = Trace(T L GT R G^, where T L , R = i(E L , R -Y,l jR ) 
are the contact spectral functions. 

In this study C 2 4-Hi2, C 54 Hi Sl C 96 H 2 4 and Ci5o-ff 3 o 
complexes have been considered in their pure, defected 
(with a single vacancy) and hydrogen functionalized 
form. All structures have been relaxed by DFT molecular 
dynamics while relaxation information is also used by the 
EH method. In the case of TB an ideal reconstruction 
of the molecular structure is considered since the latter 
does not account for interatomic distances. Finally, for 
the quantum transport calculations the islands are placed 
within two semi-infinite Ait(lll) metallic planes (directly 
considered in the case of EH, appropriately fitted in the 
case of TB[33]) in a molecular bridge configuration. 



III. COMPARATIVE ANALYSIS OF THE 
ELECTRONIC STRUCTURE BETWEEN 
FIRST-PRINCIPLES AND SEMIEMPIRICAL 
METHODS 

A proper treatment of quantum transport modeling 
has to take care of both quantitative and qualitative as- 
pects of the electronic structure of a molecular system. 
In this sense, if the value of the HOMO-LUMO gap is a 
quantitative feature, the form of the HOMO and LUMO 
wavcfunctions, or similarly, the local density of states 
of the structure for energies near the HOMO/LUMO 
states are qualitative characteristics. It can be ar- 
gued that in terms of conduction, although the former 
can influence scaling, the latter can affect the shape of 
the current-voltage curve. In addition, higher-bias con- 
duction requests accuracy for entire conduction/ valence 
bands. Such considerations imply that a proper descrip- 
tion of the local density of states of a molecular system 
by means of a quantum chemical method can be fun- 
damental for the correct modeling in terms of quantum 
transport. This section examines electronic configuration 
aspects one by one. 



A. Geometry relaxation 

From a numerical point of view, distance between the 
atomic sites influences the electronic structure of a molec- 
ular system by affecting both overlap and Hamiltonian 
matrix elements. Geometry relaxation with the DFT 
method has shown that near the island edges, complex 
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FIG. 2: HOMO-LUMO gaps of pure (a) and defected (b) 
Cg n 2_H6n quantum dots by means of DFT, EH and TB. 

distance polymerization effects can be observed that tend 
to periodically increase/decrease C-C bonding for about 
5% from the equilibrium distance (C — C equ u « 1.42A). 
Such feature tends to propagate also in the inner parts 
of the clusters albeit in a continuously decreasing extent, 
while only the central benzene rings result having equal 
interatomic C distances. Distance shortenings due to hy- 
drogen passivation have been also observed in the case 
of graphene nanoribbons[ll|, whereas the effect there is 
localized near the edges. From a methodological point of 
view, variable bond-lengths can have a practical conse- 
quence in the parameterization of the hopping integral for 
the TB method, where for accuracy's sake an evaluation 
of each atomic pair distance should take place prior to the 
assignment of the integral value (in this sense the method 
becomes similar to the single 7r-orbital Huckel one). If 
such information is not available inaccuracies in the TB 
Hamiltonian can occur. It can be argued that futhcr 
complications in the correct estimation of interatomic 
distances have to be considered in the case of interaction 
between the molecular structures and a substrate (e.g. 
for the reproduction of laboratory conditions |38|-|40|). In 
this case combined melecular/substrate atomistic mod- 
elling should enlight interface bonding interactions for 
both the chemical and the structural characteristics. 



B. Energy levels 

A well-known aspect of geometrical symmetry is the 
presence of orbital degeneracies. Such feature is cap- 
tured by all methods for the pure structures, from DFT 
to TB. Eventually, symmetry breaking events (e.g. the 
presence of a single vacancy) lift such degeneracies and 
split the respective energy levels. In this subsection we 
visualize characteristics of the eigenvalue spectrum in a 
full quantum scale, in terms of energy gaps and state 
alignment over the energy axis. Figure [2] shows HOMO- 
LUMO gaps for pure and defected structures by means 
of DFT (both 3-21g and STO-3G), EH2-spd, EH-sp and 
TB, while EH2-sp results are similar to EH2-spd and are 
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not shown. Regarding the pure structures the following 
observations can be made: (i) the gap value given by the 
DFT is bigger than those obtained by the semiempiri- 
cal methods. This aspect is not directly related to the 
exchange-correlation absence in the semiempirical cases 
(since parameterization can take place on the basis of ab 
initio calculations) but with the chemical environment 
considered for the parameterization, which usually con- 
siders bulk structures where the effect of confinement 
cannot be evaluated. It is moreover interesting to ev- 
idence that the minimal basis set in the DFT method 
slightly overestimates this value with respect to the 3-2 lg 
case, (ii) The EH2-spd method with its orbital founda- 
tion and sp 2 -hybridized calibration approximates better 
the DFT results with respect to the other semiempirical 
methods, (iii) Albeit the clear difference in terms of the 
methodology, results by EH-sp and TB are very similar. 
It should be noted here that a rough method to bring 
semiempirical models closer to first principles calcula- 
tions is by globally fitting the hopping integral (or simi- 
larly the Wolfsberg-Helmholtz K constant for EH). How- 
ever in this case the parameterization looses any meaning 
outside the designated geometrical environment. 

The validity of the results obtained by the first- 
principles B3LYP/32-lg model for these structures with 
respect to the semiempirical methods is tested by direct 
comparison with experimental data on the fundamental 
optical gap of the Ci±H\2 island. For this purpose we 
follow a TD-DFT approach for the calculation of the ex- 
citation energies of this cluster. The calculated value 
for the fundamental optical gap is E opt = 3.18eV, which 
is in a good agreement with the ex per imental value of 
E op t — 3.29eU measured in Ref. 38] [49]. The difference 
between the ground state (Ehomo-lumo = 4.13eV) 
and the excited state (E opt — 3.18eU) is also consistent 
with experimental measurements for similar structures, 
since the exciton binding energy in these complexes gives 
rise to a 0.5 — leV reduction of the optical gap with re- 
spect to the ground state HOMO-LUMO gap[4l|. 

Moving on to the defected structures (with a single va- 
cancy in the central benzene ring) an expected reduction 
of the gap value can be observed, while DFT basis set 
differences become more pronounced. For the semiem- 
pirical methods the picture changes qualitatively only 
in the EH2-spd case. Considering spin-degeneracy, the 
EH2-spd parameterization assigns the HOMO and the 
LUMO states to two quasi-degenerate levels prior to the 
real gap, which in this case is represented by the LUMO 
and LUMO+1. A detailed study of the corresponding 
eigenvectors (with respect to eigenvectors given by the 
DFT and EH-sp) shows that in the case of the defected 
structures an incorrect state is inserted at the energy 
axis inside the energy gap. In this sense, although the 
EH2-spd/EH2-sp parameterizations demonstrate overall 
optimal characteristics (e.g. see ref. 35] for a study on 
carbon nanotube band structure and the next subsection 
for pure clusters), a careful use might be necessary for 
transport calculations in defected coronene systems. 




-1 o 1-1 1 



E - E F (eV) 

FIG. 3: Density of states around the Fermi level for pure 
(black line) and defected (red line) n — 5 islands by means of 
(a) DFT (STO-3G), (b) DFT (3-21g) (c) EH-sp and (d) TB. 
In all cases, a small smearing has been applied. 

Another important aspect that regards energy eigen- 
states is their alignment over the energy axis, moreover 
when local alterations of the symmetry 'break' the ideal 
atomic structure. For few-atom molecular complexes a 
neat way to visualize this is with their density of states 
spectrum as a function of energy. Figure [3] plots DOS 
functions for pure and defected n = 5 islands by means 
of DFT, EH-sp and TB. In this case the lifting of the 
symmetry-induced degeneracy inserts states that tend to 
'shrink' the band gap. In the case of TB, a state always 
appears in the center of the pure structure's gap. This 
well-known effect has its origin at the bipartite nature of 
the honeycomb lattice, where the presence of a single va- 
cancy in one of the two sub-lattices inserts a zero energy 
mode at the Fermi level of the system (i.e. at energy 
E=0)[42j]. The key issue though arising from figure is 
that for the more sophisticated methods the HOMO state 
is not located in the center of the pure structure's gap, 
but shifted towards lower energies next to the valence 
band. Such feature is captured by both DFT and EH-sp, 
although in a quantitative disagreement. The analysis 
therefore implies that TB gives a rigid picture of the gap 
state with respect to more sophisticated models. Under 
this perspective, this study will try to affront the problem 
by introducing a further parameterization for the point 
defect (see section \FV\i . 

C. Qualitative evaluation: molecular orbitals 

A most important aspect for transport in nanostruc- 
tures is the availability of states (cither full or empty) 
within the conduction window. In conjunction, a very 
important factor for the correct treatment of conduction 
are the eigenvectors that correspond to these states, from 
which topological features can be deduced (e.g. localiza- 
tion, polarization etc.). In this context the semiempirical 
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FIG. 4: Molecular orbitals for the HOMO-1, HOMO, LUMO 
and LUMO+1 states of the n — 2 complex by means of DFT 
(upper) and EH-sp (lower). 



FIG. 5: Molecular orbitals for the HOMO-5, HOMO-4, 
HOMO-3 and HOMO-2 states of the n — 2 complex by means 
of DFT (upper) and EH2-sp (lower). 



methods have been evaluated on the basis of DFT results 
for pure, defected and hydrogen functionalized struc- 
tures. Results shown here are for the DFT 3-21g, EH-sp, 
EH2-sp and TB models. A qualitative correspondence 
has been obtained for the DFT STO-3G and the EH2- 
spd bases with respect to their method counterparts [50] . 



1. Pure structures 

TB: In the absence of real atomic orbitals and with 
the restriction of its limited basis set the TB method 
demonstrates a progressive wavefunction descriptive ca- 
pacity from smaller to larger complexes. In particular 
TB eigenvectors for conduction/ valence eigenstates have 
a poor resemblance with the respective DFT ones for 
the C24H12 molecule, while similarity becomes gradually 
better for bigger structures. This behavior is due to the 
C-C bonding distance polymerization features discussed 
earlier, which have a higher impact for the smaller com- 
plexes. Arriving at the n = 5 complex, DFT-TB match- 
ing becomes adequate, hence, the TB method is qualified 
for the electronic structure description of these structures 
with a respective number of C atoms and onwards. 

EH-sp: The method is in a qualitative agreement 
with the DFT one only for the energy degenerate 
HOMO/HOMO-1 and LUMO/LUMO+1 pairs for all 
studied structures (see fig. 0] for the coronene molecule). 
Moving away from these states towards the valence band 
accuracy is lost, not in the form of the wavefunctions 
whereas in the correct order that these appear. Conduc- 
tion band description results poor. 

EH2-sp: Matching between EH2-sp and DFT wave- 
functions is excellent for all pure islands and for both 
valence and conduction energy zones (e.g. see figure [5] 
for the valence band of C^Hx-i)- It is evident here that 
the chemical environment in which the parameterization 
has taken place (bulk graphene) and the double-exponent 
Slater orbitals play a crucial role in the representation of 
correct molecular orbitals. It is also interesting to note 
the mismatch in the results obtained by the EH2-sp and 
the EH-sp parameterizations, even if the quantum chem- 
ical method is the same. 



Finally, magnetism issues that could arise due to the 
presence of zigzag terminated edges in these complexes 
are not confirmed by the the ab initio calculations, con- 
trary to zigzag GNRsjiTl], [l2|, |43[ or coronene islands 
above a critical size (based on mean-field Hubbard model 
calculations (3). It can be therefore stated that the 
absence of a self-consistent exchange evaluation in the 
semiempirical methods does not compromise the ob- 
tained results with respect to the DFT case in the present 
study. 



2. Defected structures 

The presence of a single vacancy in these islands pro- 
vokes a distortion of the atomic structure since the re- 
maining a dangling bonds tend to recombine by leaving 
their equilibrium positions. On the other hand, cr-orbital 
energies are too far away from the HOMO-LUMO states 
and do not contribute to the formulation of the respective 
wavefunctions. Most importantly, apart from a symme- 
try breaking effect in topological terms, the presence of 
the vacancy imposes a localization of the wavefunctions 
that correspond to the various eigenstates, making such 
complexes 'sensitive' to the positioning of a nanoprobe. 
In terms of the various methodologies we have obtained: 

TB: The lack of information concerning distance in 
the TB method is even more important for the defected 
structures, where the smaller the structure the higher 
is the effect of the vacancy on its deformation. In this 
sense the TB method with its standard parameteriza- 
tion is inadequate for the description of the electronic 
structure of these complexes, whereas like in the case 
of pure structures, description gradually betters as the 
complexes grow, with the following particularities: (i) 
HOMO wavevectors present a succession of zero and non- 
zero values for neighboring atomic sites and (ii) for even 
n-indexed molecules (n — 3, n = 5 etc.) the hexagonal 
edge that corresponds to the defected site presents atoms 
with zero LDOS for the HOMO eigenstate (fig. [B}. This 
last observation is crucial in terms of transport modeling 
and its implications will be discussed in the next section. 

EH-sp: It describes better the valence (HOMO, 
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FIG. 6: Highest occupied molecular orbitals for a) n—2, b) 
n=3, c) n=4 and d) n=5 islands with a single vacancy by 
means of DFT (left), EH-sp (middle) and TB (right). The TB 
orbital representation is purely demonstrative by assigning 
Slater type p z orbitals with the same parameters as the EH 
method. 

H0M0-1, H0M0-2) than the conduction band. More- 
over, moving towards the bigger structures accuracy is 
increased and for the n — 5 island description becomes 
adequate for the valence and discrete for the conduction 
band. 

EH2-sp: The main drawback of this parameterization 
has to do with the presence of an incorrect HOMO eigen- 
state as discussed in the previous subsection. Overall 
it offers a valid alternative to the DFT results, on the 
other hand though, the importance of the HOMO state 
in terms of conduction modeling requires attention in its 
use in defected graphene environments. 

Finally a remark on the C3 point symmetry of the 
HOMO wavefunction around the vacancy should be 
made[42j (feature that is captured by all methods), where 
clearly a non-zero magnetic moment arises (also ob- 
tained with Hartree-Fock-based calculations on the same 
complexes (18|)- 



3. Hydrogen-functionahzed structures 

Results for the complexes where the defected site has 
been functionalized by a hydrogen atom that saturates 
one a dangling bond do not differ substantially from 
their nonfunctionalized counterparts. Here the role of 
hydrogen slightly influences the structure's geometrical 



relaxation whereas wavefunctions are similar to non- 
passivated molecules with defects, as er-orbital energies 
are too far away from the zone of interest for conduc- 
tion. Consequently the discussion made in the previous 
subsection is valid also in this case for the methods that 
directly account for the presence of hydrogen (DFT, EH). 



D. Conclusion 

Semiempirical models in graphene-based quantum dot 
structures can be successfully used within a certain 
framework that is established by their quantum chem- 
ical limitations. The Extended Hiickel method with its 
real-orbital foundation can cope with a great number of 
qualitative features, whereas the role of parameterization 
proves to be fundamental. In this sense EH2-spd/EH2- 
sp are excellent alternatives to DFT for pure coronene 
structures, whereas defected/functionalized complexes 
are more appropriately treated by the EH-sp model. On 
the other hand, TB with the standard parameterization 
can be used for the study of conduction in large defect- 
free systems, while modeling remains a challenge for de- 
fected complexes since results appear too 'radical'. In 
this sense a further parameterization of the defected site 
within its particular topological environment is necessary 
for the correct estimation of the electronic structure, ar- 
gument that will be treated in the next section. 



IV. QUANTUM TRANSPORT 

The study's bottom line is to efficiently model trans- 
port phenomena on graphene-based quantum dot sys- 
tems respecting the chemical aspects that arise due to 
the particularity of the chemical/geometrical environ- 
ment. The importance of a proper description of the 
electronic structure on conduction can be better appre- 
ciated in the case of defected islands, where the pres- 
ence of the vacancy is a reason for topological asym- 
metries also on the formation of the molecular orbitals 
(see fig. IS]). In this section, a numerical analysis takes 
place for the defected n = 5 complex initially with the 
EH-sp method, while results are used for a critical eval- 
uation of similar calculations made with the TB model. 
Two equivalent molecular bridge configurations are used, 
where the source-device-drain geometry differs only in 
the position of the contacts with respect to the vacancy 
site (fig. [7]) . In detail, two opposite edge corners of the 
aforementioned dot have been inserted between two semi- 
infinite Au(lll) metallic planes, which model the metal- 
lic probes of an atomic force microscope. The contact 
Hamiltonian is also written within the EH theory using 
an appropriate spd basis 33]. A prerequisite of equiva- 
lence for the contact bonding between the two configura- 
tions is explicitly requested for an evaluation of trans por t 
without geometrical or bond strength implications [33] . 
Here, contacts are 1.7A distant from the edge C atoms 
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(avoiding strong invasiveness) and are geometrically sym- 
metrical with respect to the molecular structure. The 
transmission probabilities obtained for the two configu- 
rations have a distinct character, whereas differences are 
not fundamental for the conduction characterization of 
the system. Namely, both configurations give a non-zero 
transmission value corresponding to the HOMO state, 
whereas differences exist both in the valence and con- 
duction band. The divergences can only be attributed to 
the different geometrical positions of the contacts with 
respect to vacancies that reflect unequal interface chemi- 
cal bondings due to orbital localization phenomena. The 
minor impact of such phenomena on the conduction char- 
acteristics is driven by the real atomic orbital founda- 
tion of both contact and device wavefunctions that con- 
stitute bonding interactions that exceed next neighbor 
distances. Therefore, e.g. if local disorder provokes a 
nullification of the LDOS at the contact-device interface 
at a certain energy, transmission is still possible if this 
zero LDOS expands in a smaller area than that of orbital 
overlap between contact and further device atoms with a 
finite LDOS. The same concept can be described from a 
quantum mechanical perspective, where the presence of 
the contacts induces a constant perturbation on the bare 
device's Hamiltonian and the effective Hamiltonian now 
writes: 

Heff = H o + Hql + H r (2) 

Here Hq corresponds to the molecular Hamiltonian in 
the absence the contacts and i?oL.O-R ar e the Hamiltonian 
components that arise due to the interaction between the 
device and the left/right contact. According to EH the- 
ory for the localized HOMO state we get: 

< ^homo\H l\^homo (3) 

and 

< ^homo\H r\^ HO mo (4) 

This finite value of both integrals makes transport plau- 
sible from the HOMO state. 

Albeit its phenomenological simplicity, the TB descrip- 
tion of conduction in the same structures as before gen- 
erates complexity in the interpretation of the obtained 
results. The critical points are two: (i) the first has to 
do with the zero mode introduced by the vacancy at the 
Fermi energy level. The corresponding HOMO wavefunc- 
tion, and equally the LDOS at E = 0, have a succession of 
finite and zero values for next-neighbor atoms, that is, for 
each C atom with a finite LDOS value the three nearest 
neighbor atoms have a zero value and vice versa. More- 
over the n = 5 structure (like all even n-indexed ones) has 
a hexagonal side with LDOS = 0, as discussed in the pre- 
vious section. Therefore, for the standard-parameterized 
first-neighbor TB model this state represents the respec- 
tive molecular orbital in a rigid way, contrary to EH and 
DFT. (ii) The second issue reflects TB interface bonding 
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FIG. 7: Transmission as a function of energy by means of EH- 
sp for the defected n — 5 complex, for two equivalent contact 
configurations that differ only in the position with respect to 
the defected site. 




Energy (eV) 



FIG. 8: Density of states and Transmission probability of the 
n = 5 complex by means of the TB method for nonparam- 
eterized (a) and parameterized (b) vacancy values (Eon — 
site vac = lOeV, to = 1.9eV). The DOS figures are repre- 
sented with a small Gaussian smearing. 

issues between a device and the metallic leads that in the 
next-neighbor context present a strongly-localized char- 
acter (e.g. only two C atoms in our case are allowed to 
chemically interact with the leads). In this case, if the 
metallic contacts form bonding interactions exclusively 
with zero LDOS carbon atoms (configuration 2 in our 
case), the HOMO eigenstate will not contribute to the 
conductivity of the system, yielding a zero transmission 
probability at that energy. Indeed, figure [5] shows trans- 
mission as a function of energy for the n — 5 complex for 
the two contact configurations presented before, where a 
finite transmission probability for the HOMO state ap- 
pears only in the first case, whereas clearly no conduction 
takes place through this state for the second. In terms 
of expectation values, for the HOMO eigenstate of the 
second configuration we now get: 

< ^ homo\H r\^ homo >= (5) 

This blocked conduction channel reflects the extreme 
manifestation of wavefunction localization obtained by 
TB and comes to contrast with EH results. It is there- 
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FIG. 9: Schematic LDOS representation of the n = 5 com- 
plex by means of the TB model for two different parame- 
terizations of the vacancy site, a) to = OeV within the va- 
cancy and the neighboring sites, and b) to = 1.9eV and 
Eon — site vac ~ WeV. The radius of each circle is propor- 
tional to the amplitude of the LDOS value on that atomic 
site. 

fore fundamental that a realistic modeling has to take 
into account that the tails of the contact wavefunctions 
penetrate the body of a molecular device for several A 
before they decay. 

A possible way to affront the aforementioned problems 
is by introducing a further parameterization of the va- 
cancy site. Figure M shows a schematic LDOS real-space 
representation of the n — 5 island by means of the TB 
method for a non-parameterized and a parameterized va- 
cancy site. The vacancy parameterization takes place by 
assigning a finite E = fOeU energy on the site and a 
to = l.9eV hopping integral within this and the neighbor- 
ing atoms. The principal differentiations obtained are: 
(i) The sites where zero LDOS values corresponded for 
the HOMO level now obtain a finite, albeit small den- 
sity (not visible in fig. [SJ. (ii) The hexagonal edge 
that corresponds to the vacant site (which for even n- 
indexed molecules had zero HOMO-wavcfunction com- 
ponents) obtains also a finite LDOS value that is more 
similar to the electronic structure by means of the DFT 
and the EH methods (see figure [6]). (hi) The colloca- 
tion of the HOMO eigenstate on the energy axis is in 
E < 0, i.e. it moves towards the valance band leav- 
ing the midgap position (see figure [3]). Also in this case 
the DOS spectrum comes closer to the ones obtained by 
DFT and EH (figure [3J . Finally, changes obtained for 
wavefunctions that correspond to other than the HOMO 
level do not present particular differences from their non 
parameterized counterparts. Overall, the parameteriza- 
tion of the vacancy site permits for a clear improvement 
of the qualitative aspects of the electronic structure for 
defected molecules, respecting the chemical equilibriums 
and approaching results obtained by more sophisticated 
methods. Indeed, in terms of transport, the HOMO 
eigenstate now contributes to conduction for both con- 
tact configurations whereas a less 'radical' representation 
of the transmission probability is sketched. Apart from 
to the qualitative gains of TB calibration presented here 
though, a key conceptual issue arises. Now the vacancy 



site becomes similar to a (nominally p-type) impurity, 
since the local point potential lowers (from oo to big fi- 
nite) and the hopping integral raises (from OeV to finite) 
[42| . This consideration can have an impact on the way 
vacancies are seen in graphene-based systems, both from 
an applicative as well as from a methodological point of 
view. As a conclusion, it should be strongly stated that 
the common perception of treating vacancies in graphene- 
based systems (zero-energy modes) is not confirmed in 
this study, whereas an impurity-like behavior has been 
obtained. 



V. DISCUSSION 

One important aspect of the understanding of impu- 
rity induced disorder in graphene is the possibility of a 
controllable band gap tailoring for semiconductor appli- 
cations [43, Eg] ■ This study evidences that vacancies in 
graphene complexes actually behave as impurities. Such 
consideration can have a big practical impact on the engi- 
neering of mobility gaps in graphene-based systems since 
vacancies are easier to obtain (e.g. by ion irradiation [471 ]) 
than actual p or n-type doping. Here we have attempted 
to affront modeling issues for pure and defected graphene 
quantum dot islands keeping in mind that the desired 
computational efficiency for the simulation of large sys- 
tems should not be in contrast with chemical accuracy. In 
this sense a multiscale approach has been introduced with 
the scope to identify merits and limitations of semiem- 
pirical approaches within a designated chemical environ- 
ment prior to their use for the calculation of quantum 
transport. Model confrontations have demonstrated that 
no perfect matching exists between the results obtained 
by the ab initio on the one hand and the semiempirical 
approaches on the other. The extended Hiickel method 
with its real-orbital foundation manages to capture a 
wide set of qualitative aspects of the systems, which qual- 
ify it as an appropriate method for quantum transport 
calculations in graphene-based environments. Moving to- 
wards computational efficiency, the tight-binding model 
has confirmed its authoritativeness for pure large-scale 
structures, whereas when structural defects have to be 
accounted a further parameterization of these sites needs 
to be considered. Unlikely, such tuning cannot be generic 
for all types of complexes/defect-types since chemical en- 
vironment influence can be fundamental. E.g. we have to 
note that by only adding a second vacancy in the vicin- 
ity of the same triangular sublattice of the honeycomb 
structure, hybridization between the two modes can take 
place. In this sense, a model evaluation of TB by a more 
sophisticated method should ideally take place prior to 
its use in disordered graphene-based systems. It should 
be finally pointed out that both theoretical and experi- 
mental attention should be paid to strongly defected sys- 
tems where topological disorder can be a reason for wave- 
function localizations, whose influence on the electronic 
properties can be important. 
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We have confirmed the good TD-DFT estimation of the 
optical gap with respect to the experimental value also 
by calculating excited state energies within the Config- 
uration Interaction Singles (CIS) Method. The optical 
gap given by the latter is E op t = 4.19eV , which is by far 
bigger than the experimental and the TD-DFT value. 
The surprisingly accurate correspondence in terms of 
molecular orbitals obtained for the 3-21g and STO-3G 
bases within the DFT scheme indicates that the mini- 
mal basis set can be used in disordered graphene-based 
systems with only small quantitative compromises (see 
paragraph IIHI B) . 



